Incorporating the molecular gas phase in galaxy-size numerical simulations: 

first applications in dwarf galaxies 
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P^. ABSTRACT 
O 

' We present models of the coupled evolution of the gaseous and stellar content of 

■ galaxies using a hybrid N-body/hydrodynamics code, a Jeans mass criterion for the 

onset of star formation from gas, while incorporating for the first time the formation of 
^ , H2 out of HI gas as part of such a model. We do so by formulating a subgrid model 

H ' for gas clouds that uses well-known cloud scaling relations and solves for the HI<->H2 

balance set by the H2 formation on dust grains and its FUV-induced photodissocia- 
tion by the temporally and spatially varying interstellar radiation field. This allows 
the seamless tracking of the evolution of the H2 gas phase, its precursor Cold Neutral 
Medium (CNM) HI gas, simultaneously with the star formation. An important advan- 
tage of incorporating the molecular gas phase in numerical studies of galaxies is that 
the set of observational constraints becomes enlarged by the widespread availability of 
H2 maps (via its tracer molecule CO). We then apply our model to the description of 
the evolution of the gaseous and stellar content of a typical dwarf galaxy. Apart from 
their importance in galaxy evolution, their small size allows our simulations to track 
the thermal and dynamic evolution of gas as dense as n ~ 100 cm~^ and as cold as 
Tk ~ 40 K, where most of the HI — > H2 transition takes place. We are thus able to 
identify the H2-rich regions of the interstellar medium and explore their relation to the 
ongoing star formation. Our most important findings are: a) a significant dependence 
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of the HI H2 transition and the resultant H2 gas mass on the ambient metallicity and 

the H2 formation rate, b) the important influence of the characteristic star formation 
timescale (regulating the ambient FUV radiation field) on the equilibrium H2 gas mass 
and c) the possibility of a diffuse H2 gas phase existing well beyond the star-forming 
sites where the radiation field is low. We expect these results to be valid in other types 
of galaxies for which the dense and cool HI precursor and the resulting H2 gas phases are 
currently inaccessible by high resolution numerical studies (e.g. large spirals). Finally 
we implement and briefly explore a novel approach of using the ambient H2 gas mass 
fraction as a criterion for the onset of star formation in such numerical studies. 

Subject headings: galaxies: numerical simulations -galaxies: dwarf irregulars - galaxies: 
star formation - ISM: molecular gas - atomic gas 

1. Introduction 

Stars form in molecular clouds, and most stars form in the large complexes of Giant Molecular 
Clouds (GMCs) . A general theory to predict the location and amount of star formation in galaxies 
does not exist yet. The global drivers of star formation are generally sought in large scale processes 
that can form concentrations of HI gas (Elmegreen 2002), for example gravitational instabilities in 
disk galaxies, or compression by galaxy interactions or mergers. These concentrations of neutral 
gas are then assumed to be the sites of molecular cloud and ultimately star formation. Numerical 
simulations of galaxies have followed this lead and base their star formation model on the local 
density directly using a Schmidt law (Mihos & Hcrnquist 1994; Springel 2000), or gravitational 
instability and a Jeans mass criterion (Katz 1992; Gcrritsen & Icke 1997; Bottema 2003). To date 
none of these efforts has incorporated the emergence and evolution of the H2 gas along with its 
precursor Cold Neutral Medium (CNM) HI phase, except in a semi-empirical and somewhat adhoc 
fashion (Semelin &; Combes 2002). The need to do so in a more self-consistent manner is now 
recognized from both numerical (Bottema 2003), and observational studies (Wong & Blitz 2002). 

The interstellar medium (ISM) consists of gas of wide ranging properties, from cold, dense 
molecular clouds to the cold and warm HI medium (CNM, WNM), as well as the hot ionized 
medium intermixed in fractal-like structures. The physical processes governing the state of the 
ISM have been progressively identified in the last few decades, but their precise working is still an 
area of active research (Vazquez-Scmadcni 2002). Most galactic-scale simulations include only one 
phase, the WNM HI (Katz 1992; Navarro & White 1993; Springel 2000). Some authors include 
more physics and consider a two-phase medium (Gerritsen &; Icke 1997; Gerritsen Sz de Blok 1999). 
Semi-empirical models of a multiphase ISM (Springel &; Hernquist 2002; Andersen &; Burkert 2000) 
suffer from the inclusion of many poorly understood parameters needed to describe the interaction 
of the various phases, which also forces such models to include serious simplifications (e.g. a con- 
stant and uniform heating of the gas) in order to keep the problem within the current computing 
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capabilities (Semelin & Combes 2002). Thus the robustness of such semi-empirical approaches 
in modeUng real galaxies is rather limited. Moreover, apart from the obvious fallacy of forming 
stars out of atomic gas, the absence of molecular gas in the simulations inhibits the comparison 
of the model galaxies with an extensive body of observational data, namely the distribution of 
H2 in galaxies (as deduced via its tracer molecule CO). The latter is observationally well-studied 
(Engargiola et al. 2003; Mizuno et al. 2001; Heifer et al. 2001; Regan et al. 2001), it can offer ad- 
ditional constraints on the model galaxies, and the same is true for its observed empirical relations 
with other galactic components e.g. the young stars. For example the Schmidt law that links local 
gas content to the star formation rate has been demonstrated to be a much tighter relation for 
molecular, rather than HI gas (Wong & Blitz 2002). 

Some pioneering work that includes H2 in spiral galaxy simulations has been done by Hidaka &: 
Sofue (2002) but they did not consider a time-dependent HI — > H2 transition, which as we will argue 
in this work, is essential. Our time-dependent treatment of the HI H2 transition, implemented 
here in an N-body/SPH code, can be incorporated into other types of code that include the ISM 
evolution thus enabling the resulting galaxy simulations to track the gas phase most relevant to star 
formation, the molecular gas. Our work fills two important "gaps" in the published literature, 
namely: a) it makes a connection between galaxy-sized simulations and molecular cloud theory 
(which lies behind the emergence of the power laws observed in the cool ISM), and b) connects large 
scale instabilities and the appearance of CMC-type complexes. The latter allows the integrated 
study of star formation and the H2 gas in the dynamic setting of realistic galaxy models that include 
effects of e.g. spiral density waves, self-propagating star formation, galaxy interactions and mergers, 
and a temporally and spatially varying interstellar radiation field. Our presentation proceeds from 
a brief exposition of the H2 formation/destruction theory and the cloud power-laws underlying 
our sub-grid physics assumptions, on to its implementation in our numerical simulations while 
discussing the relevant details of the code we use. Finally we present first results from simulations 
of typical dwarf galaxies, and conclude by outlining future work using the H2-tracking numerical 
models. 



2. H2 formation and destruction 

The formation of H2 in galaxies has already been studied by Elmegreen (1989, 1993), where 
the important role of ambient metallicity and pressure has been described. Subsequent efforts to 
incorporate his approach into analytical models of galactic disks (Honma et al. 1995) or numerical 
simulations (Hidaka & Sofue 2002) have adopted stationary models, namely once an H2 formation 
criterion is satisfied the HI H2 transition is set to be instantaneous. The time-dependence of 
several factors affecting this transition (e.g. the ambient H2-dissociating FUV radiation field) and 
the various gas heating and cooling processes were not considered, seriously restricting the ability 
of such models to track the evolution of the ISM and particularly its H2 gas phase in a realistic 
manner. This is because the conditions affecting the HI<-^H2 equilibrium can vary over timescales 
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comparable or shorter than that needed for the latter equihbrium to be reached. 

Dust grains, when present, are the sites where H2 forms. Its formation rate Rf for gas with 
temperature and metaUicity Z can then be expressed as 

1 /2 

Rf = l<ra{vn)Snln, = 3.5 x IQ-^Vz (x^) ShTh^ cm^ s-\ (1) 

(e.g. Hollenbach, Werner & Salpeter 1971; Cazaux Sz Tielens 2002, 2004). This equation expresses 
the H2 formation rate as the rate at which HI atoms with mean velocities (vh) collide with inter- 
stellar dust particles with effective surface a^, multiplied by the probability SH(Tk) they stick on 
the grain surface and the probability 7H2 they eventually form an H2 molecule that detaches itself 
from the grain. The grain surface (and thus Rf) is assumed to scale linearly with metallicity 
Z, namely = dgng/n = 4.9 x lO^^^Z cm^ (with (jg the grain geometric cross section and with 
the ratio of grain density to total hydrogen density Ug/n oc Z). Laboratory experiments usually 
constrain only the Sh7H2 product for various temperature domains rather than provide information 
for each function separately (e.g. Pirronello et al. 1997; Katz et al. 1999). The theoretical study 
of Buch & Zhang (1991) yields a function Sh = [1 + (keTk/Eo)]^^ (with a characteristic energy 
scale Eo/ke ~ 100 K obtained from fitting results from molecular dynamics simulations) valid for 
Tk < 300 K, which we adopt in the present work (but see also Cazaux &; Spaans 2004 for the most 
recent views). 

We incorporate all the uncertainties of Rf into a single parameter /x. These uncertainties stem 
mainly from the adopted probability functions and the effective surface for H2 formation. For 
example, the effective H2 formation surface may be bigger than implied by the visual extinction 
cross-section a^. The canonical formation rate Rf = 3 x 10-^''' cm^ s-^ (Jura 1974, 1975) for typical 
CNM HI gas (Tk ~ 100 K), corresponds to = 3.5 but values of ~ 17 — 18 are not excluded. 
Early suggestions for such high /j, values emerged from the study of rovibrational infrared lines of 
H2 in reflection nebulae (Sternberg 1988), and observations of intense H2 rotational lines with the 
Infrared Space Observatory (ISO) in photodissociation regions (Habart et al. 2000; Li et al. 2002). 

The timescale associated with H2 formation is then (e.g. Goldshmidt Sz Sternberg 1995) 

(rp \ —1/2 _-|^ 

iOolc) (10^)" bZSH(T.)l-' yrs, (2) 

where ni and Tk are the HI density and temperature. For ni ~ 50 cm"'^ and Tk ~ 100 K, typical 
of the Cold Neutral Medium HI gas out of which H2 clouds form, rf ~ 10^ yrs. This is comparable 
to the timescales of a wide variety of processes expected to fully disrupt or otherwise drastically 
alter typical molecular clouds and their ambient environments. Some of the most important ones 
are star formation, with the disruptive effects of O, B, star clusters (Bash et al. 1977), turbu- 
lent dissipation (Stone et al. 1998; Mac Low et al. 1998), and inter-cloud clump-clump collisions 
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(Blitz & Shu 1980). The mean FUV field driven by the evolution of continously forming stellar 
populations throughout a galaxy evolves over similar timescalcs (e.g. Parravano et al. 2003), and 
the same seems to be the case for the ambient pressure environment and its perturbations by pass- 
ing SN-induced shocks (e.g. Wolfire et al. 2003). Hence a realistic model of the HI—>H2 transition 
in galaxies must he time-dependent. 

With this necessity in mind, we will first discuss the equilibrium fraction of molecular gas in 
the ISM, before presenting a time dependent version of the model. 

2.1. The equilibrium molecular fraction 

The equilibrium molecular gas fraction per cloud, under a given ambient FUV field, can be es- 
timated by considering the formation/destruction equilibrium for a plane parallel cloud illuminated 
by a radiation field Gq. In equilibrium formation balances destruction, namely 

Rfnm = Gokof(N2)e-^n2, (3) 

which must hold at any depth. The densities ui and n2 denote the HI and H2 densities, while the 
total hydrogen density is n = ni + 2n2. The H2 dissociation rate ko = 4 x 10~^^ s^^ is normalized 
for an ambient FUV field in units of the Habing (Habing 1968) field value (Go = 1). The factor 
f(N2) is the normalized H2 self-shielding function which describes the decrease in dissociation rate 
due to the FUV absorption by the molecular column N2. Furthermore, r = (t(Ni -|- 2N2) is the 
FUV optical depth due to grain extinction, with a dust FUV absorption cross-section a = ^FUV^d 
(Cfuv = 2 — 4). Eq. 3 can be converted to a separable differential equation for Ni and N2 (e.g. 
Goldshmidt &: Sternberg 1995) which, when integrated by parts for a uniform cloud, yields a total 
HI column density of 

This is the total steady-state HI column density resulting from the FUV-induced H2 photodissocia- 
tion of a one-side illuminated uniform cloud, and can be taken as a measure of the total 'thickness' 
of the HI layer. The factor $ is an integral of the self-shielding function over the H2 column, which 
encompasses the details of H2 self-shielding. For f(N2) ~ (N2/Nch)~'^ it is 

I'QO 

$ = 2^^-^ (Ncha)'^ / X- V^dx = 6.6 x lO-60FzV2 ^i/J^. (5) 
Jo 

Its numerical value was obtained for k = 1/2 and a characteristic column density of Nch = 1-75 x 
10^1 cm-2 (Jura 1974). 
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The equation of the HI ^ H2 balance still yields an analytical expression for Ntr(HI) for 
clouds with density profile: n(r) = nc(r/R)~^ (R: the cloud radius). Such a profile results from 
cloud models using a logotropic equation of state (McLaughlin & Pudritz 1996), hence we will 
refer to them as "logotropic clouds," but the reason for considering such profiles here is to obtain 
a rough representation of sub-resolution, non-uniform, GMC density structure. The corresponding 
expressions for the equilibrium transition column density etc. for logotropic clouds are given in the 
Appendix. 

The visual extinction corresponding to a transition column Ntr(HI) is given by (using also the 
expressions for Rf , $ from Eqs 1 and 5) 



Af ) = 1.086Cp^v In 

For a spherical constant density cloud the H2 gas mass fraction contained "below" the HI transition 

ftr) 

zone marked by Ay is then given by 



1 -L Go / ^FUv\ 

MSH(Tk) V ZTk J 



1/2 



( V 

V135 cm-3y 



(6) 



^ M(H2) _ / 4Af) A 
- Me - y 3(A.) ) ' 

where (Av) is the area-averaged extinction of the cloud measured by an outside observer (fm = 
for A^''^ > 3/4(Av)) . 

We consider fj^ to be a measure of the local molecular gas content of the ISM, the implicit 
assumption being that most of the mass of the cool HI/H2 gas can be ultimately traced to such 
cloud structures. For CNM HI and the resulting H2 clouds this is a reasonable assumption and in 
certain models it is indeed predicted to be the final configuration of any initial warm (~ 10^ K) HI 
gas mass which quickly cools and fragments in an environment of constant "background" pressure 
(Chieze 1987; Chieze & Des Forets 1987). 



2.2. Sub-resolution physics: the cloud size-density relation 

From the previous discussion it becomes apparent that in order to calculate fm for a region 
of the ISM we need an expression for the typical local mean cloud extinction (A^). The crucial 
assumption we make here to derive an estimate for (Ay) is that all unresolved cloud structures where 
the HI/H2 transition takes place obey the widely observed density-size power law (Larson 1981). 
The bulk of the H2 gas is indeed observed to reside in Giant Molecular Clouds, the sites of most 
star formation activity in galaxies, and these arc observed to be virialized entities (Larson 1981; 
Elmegreen 1989; Rosolowsky ct al. 2003). Their density/size scaling relation can be derived from 
an application of the virial theorem to clouds under a given boundary pressure Pg (Elmegreen 1989) 
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and a universal (velocity dispersion)- (cloud size) relation. Then it can be shown that the average 
H density is given by 



where Rpc is the cloud radius in parsecs. The range of the normalizing constant Uq (cm~^) can 
be determined theoretically by choosing a range of plausible polytropic solutions to the cloud 
density profiles (Elmegreen 1989), or from the observed n-R relations found for galactic molecular 
clouds. The latter approach, after using the scaling relation n(H2) = ITOORp^" cm~^ reported 
by Larson (1981) (a ~ 1), yields Uq ~ 1520 cm~^. In estimating Uq from the observed n-R 
scaling relation we took into account that for the galactic midplane molecular clouds the outer 
pressure Po/ke ~ 10^ K cm~^, but the boundary pressure on the molecular part of the cloud 
Pg = Pg ~ 5Po (Elmegreen 1989). 

It is worth mentioning that the invariance of the linewidth-size relation AV(L) ~ Vo(L/pc)V2 
(which for virialized clouds yields Eq. 8) has been recently verified over an extraordinary range 
of conditions as well as for the environments within clouds (Heyer & Brunt 2004). This lends 
additional support to the use of the n-R scaling law as a non-evolving element of our assumed sub- 
grid physics. These relations are expected to break down only at very small scales where linewidths 
become dominated by thermal rather than macroscopic motions. For typical CNM temperatures 
setting AV(L) ~ AVth yields L < 0.5 pc, much smaller than the spatial resolution achievable in 
galaxy-size simulations like ours. 

The average extinction of a cloud with radius R as measured from outside is, 



where the geometric factor kg = 2/3 for a spherical cloud. After substituting (n) from Eq. 8 and 
setting Na^=i(H) = 1.85 x lO^^cm"^ we obtain 



The cloud boundary pressure is due to thermal as well as macroscopic motions (Elmegreen 
1989), and assuming the latter to be isotropic, 



i 2 

Pe = UekeTk -I- ^PeO-yel 



L.085Tk + 54 f--^^)\ 
Vkm s"^/ 



HekB, 



(11) 



where (Tyei is the 3-dimensional velocity dispersion of the gas due solely to macroscopic motions. 
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Prom Eqs 7 and 10 it is obvious that high pressure gas will tend to be molecular, a result al- 
ready known from the analytical/steady-state models of the HI — ^ H2 transition in galactic disks 
(Elmegreen 1993; Honma et al. 1995). 

We can now compute the equilibrium fm values for interstellar gas for a given (n, Tk) sampling 
the local conditions in the ISM by using Eq. 7 (or Eq. B3 for logotropic clouds). Input parameters 
are the radiation field Go, metallicity Z, velocity dispersion cxvei and our H2 formation rate parameter 
H. In Figure 1 we show the resulting equilibrium fraction fm for a number of typical values of these 
parameters. The true equilibrium molecular gas content will most likely be higher than fm since the 
real FUV field incident on GMCs is not radial and thus falls off more rapidly inside the absorbing 
cloud layers. Higher gas densities expected to exist deeper in the clouds as part of any spatial and 
density substructure/clumping in the ISM will act to raise the true fm- In that respect logotropic 
clouds are more realistic than those with uniform density. 



2.3. Time-dependence of the HI/H2 equilibrium 

The time-dependent HI/H2 transition can be approximated from the solution of Equation 22 
of Goldshmidt & Sternberg (1995), which describes the evolution of the total HI column Ntr(HI) 
in the HI H2 transition layer of a plane parallel slab of H2 gas with an impinging dissociating 
radiation field: 



^ daNUm,t) ^ ,^^^,-.N.(Hi,t) _ ^N,,(HI, t), (12) 

where Tf = l/(2nRf) (the H2 formation timescale in the fully atomic part of the cloud), and the 
dimensionless quantity r^is quantifies the local balance between H2 dissociation and formation and 
is given by 



The assumption underlying the validity of Eq. 12 is that of a sharp HI/H2 transition layer separating 
the slab into a fully atomic and a fully molecular section. In the case of steady state Eq. 12 reduces 
to 



rdise-^N-(«i)=aNtr(HI), (14) 

which does not yield Ntr(HI) as expressed in Eq. 4 since the latter does not involve the approxi- 
mation of a sharp HI/H2 transition made here. 
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In order to see at what step of deducing Eq. 4 the assumption of a sharp HI/H2 transition 
layer would yield the result in Eq. 14 we can use Eq. 2 from Goldshmidt and Sternberg (1995) 

RfndN(HI) = Gokof[N(H2)]e-2'^N(^2)e-'^^(™)dN(H2). (15) 

For a sharp HI/H2 transition zone not much HI can exist in that zone to significantly add to the 
total HI transition column density. Hence, when integrating the last equation by parts, e"'^'^^^^) 
can be replaced by e""''^*'^^^^). The latter is then treated as a constant of integration over the H2 
column density and kept out of the integral, yielding Eq. 14. We find the steady-state solutions 
provided by Eq. 14 to differ < 30% with those provided by the more accurate Eq. 4 while other 
uncertainties of our model (e.g. the effects of CNM HI gas substructure) are expected to be more 
important. The sharpness of the HI/H2 transition zone (AAy ^ 0.1 for Z=l) is due to the extreme 
self-shielding properties of the H2 molecule and has been known since the early theoretical models 
of Photodissociation Regions (PDRs) (Hollenbach, Takahashi &; Tielens 1991; Hollenbach & Tielens 
1999). 

A fully time-dependent treatment of the HI/H2 transition is provided by solving Eq. 12 and 
then using Av"^^ = ^p^y [(7Ntr(HI, t)] to find the corresponding frn(t) from Eq. 7. Although the 
approximation of a narrow transition layer is not strictly valid, the time dependent solution of the 
full integro-differential equation would be too costly for use in numerical simulations, and also the 
uncertainties inherent in our model would not justify the additional expense for the more exact 
solution. In Figure 2 we compare the solution of Eq. 12 with the full solution as calculated by 
Goldshmidt &; Sternberg (1995). The main difference is that our solution for the Ntr grows roughly 
30% slower, which gives an effect similar to an error of 30% in e.g. the grain opacity or H2 formation 
rate parameter. 

2.4. CoUisional destruction of H2 at high temperatures 

The preceding description of the HI/H2 balance is not valid for gas much warmer than its 
CNM phase where coUisional rather than FUV-induced destruction of H2 dominates. Under 
such conditions Rf ~ since once T^ ^ 1000 K, the H2 formation effectively becomes zero 
(Cazaux & Tielens 2004) and any remaining (from the FUV destruction) H2 will eventually be 
collisionally destroyed, especially at temperatures Tk > 3 x 10^. Moreover, at such high tempera- 
tures the cloud scaling laws (and hence the expressions used to estimate fm) no longer hold since 
the observed "cloud" linewidths are then almost purely thermal and thus no longer scale with size. 

In principle the coUisional destruction of a purely molecular gas will first be dominated by 
H2-H2 collision processes and later on by the more efficient H2-HI collisions. If we look at the 
coUisional destruction coefficients for these two processes we see that for the H2 collisions to be 
more important the condition 
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(16) 



must be satisfied at Tj, ~ 3 x 10^ K (Martin, Kcogh, & Mandy 1998). In the WNM gas phase such 
a high percentage of remaining H2 (remaining from the dominant FUV destruction mechanisms 
operating in the CNM and the CNM-to-WNM phases) is unhkely. We thus consider only the HI 
collisional destruction term as important and then the H2 density will be controlled by 



where 71 is the H2-HI collisional destruction coefficient. The (analytic) solution of this will 
be employed for the high temperature domain. Here we must mention that we assumed no FUV- 
induced HI/H2 spatial segregation in the WNM phase. If there is some remnant spatial segregation 
H2 would be destroyed by the much less effective H2 collisions, thus our assumption most likely 
overestimates the level of collisional H2 destruction in the WNM phase. 



We implement a module incorporating the aforementioned "recipe" to track the HI H2 
gas phase interplay in an N-body/SPH code for galaxy simulations that includes a model of the 

neutral gas phases. N-body/SPH codes are in routine use within the astrophysical community 
as tools to explore problems in galaxy evolution and formation, and their full description can be 
found elsewhere (Hcrnquist & Katz 1989; Monaghan 1992). We use the new conservative SPH 
formulation of Springel & Hernquist (2002) . The artificial viscosity used is the standard Monaghan 
(1992) viscosity. 

The code we use is a major upgrade of Gerritsen (1997) and Gerritsen & Icke (1997), mainly 
in its description of the evolution of the stellar and gaseous components. A complete description 
can be found in Pelupessy (2005), here we will highlight the aspects most relevant. 

The main feature making this code especially attractive for incorporating the HI/H2 phase 
transition is its detailed modeling of the ISM, including physics of the WNM and CNM HI phases, 
as well as the modelling of star formation and feedback. This allows us to conduct high resolution 
simulations following the ISM gas as it cools down to temperatures of T < 100 K and densities 
n ~ 100 cm~'^, conditions expected to be typical of the HI — > H2 transition phase. 

Note that the addition of our HI/H2 model does not add a major computational burden to 
the code: in practice we find that < 1% of CPU time is taken up by the molecule formation 
routines. Furthermore the necessary calculations (solving Eq. 12 using an implicit integration 
scheme) scale linearly with the number TV of particles, whereas the time consuming parts of the 



-Jl^ = -7i(Tk) (n-2n2)n2 



(17) 



3. Implementation 
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simulation (gravity and neighbour searching, which take up ^ 60% of CPU time for the simulations 
presented here) scale as N log N for an N-body/SPH code(Hernquist & Katz 1989). 

3.1. The neutral ISM model 

Our model for the ISM is similar to the equilibrium model of Wolfire et al. (1995, 2003) for the 
neutral phases and was extended from the more simplified models from Gerritsen & Icke (1997) and 
Bottema (2003). We consider a gas with arbitrary but fixed chemical abundances Xj, scaled to the 
target metallicity from solar abundances. We solve for the ionization and thermal evolution of the 
gas. The thermal evolution is solved with a predictor /corrector step as in (Hernquist &; Katz 1989). 
The various processes included in the ISM model are given in Table 1. The main differences with the 
work of Gerritsen & Ickc (1997) and Bottema (2003) arc the following: we use more accurate cooling, 
that is calculated in accordance with the chemical composition, we solve for the ionization balance 
(albeit still assuming ionization equilibrium) and we use the full photoelectric heating efficiency 
as given in Wolfire et al. (1995). Indicatively, in Figure 3 we plot the equilibrium temperature, 
ionization fraction, heating (=cooling) and pressure as a function of density. There it can be 
seen that as density varies the equilibrium state of the gas changes from a high temperature/high 
ionization state (T = 10^ K, Xc ~ 0.1) at low densities, to a low temperature/low ionization state 
(T < 100 K, Xc < 10^'^) at high densities. In between there is a density domain where the negative 
slope of the P-n relation indicates that the gas is unstable to isobaric pressure variations, the 
classic thermal instability (Field 1965). The shape of these curves and hence the exact densities of 
the thermal instability vary locally throughout the simulation, influenced by the time-varying UV 
radiation field and supernova heating. We set a constant cosmic ray ionization rate ^CR throughout 
the galaxy, assumed to be (cR = 3.6 x 10~^^ s~^. 

Although we will consider models of different mctallicitics we will not consider the effects of 
abundance gradients or enrichment here. In general the abundance gradients observed in dwarf 
galaxies are small (Pagel &; Edmunds 1981), so this is a reasonable approximation. The fact that 
we keep the metallicity constant in time means that the model as presented here is not yet suitable 
to follow the evolution over cosmological timescales or to simulate very low metallicity dwarfs (for 
the simulations presented here the evolution in Z would be less than 10% over the course of the 
simulation, even if no metals were to be lost to the intergalactic medium). 

3.2. Star formation and feedback 

The coldest and densest phase in our model is best identified with the CNM, where GMCs 
most likely form and remain embedded. We will use the simple prescription for star formation of 
Gerritsen &: Icke (1997) as our standard star formation model. It was shown to reproduce the star 
formation properties of ordinary spiral galaxies and it is based on the assumption that the star 
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formation process is governed by gravitational instability. A gaseous region is considered unstable 
to star formation if the local Jeans mass Mj, 

MJ - ? f 5^ < M,ef (18) 



6 \Gp J 

(with s the sound speed), the assumption being that structure on a mass scale Mj-cf is present in 
the ISM. Once a region is dense and cold enough that (18) is fullfilled, the rate of star formation 
rate is set to scale with the local free fall time tg: 

The delay factor fgf accounts for the fact that cloud collapse is inhibited by either small scale 
turbulence or magnetic fields (Mac Low & Klessen 2004; Shu et al. 1987). Its value is uncertain 
and we consider values fgf = 2.5 — 20. The actual implementation of star formation works as follows: 
once a gas particle is found to be unstable according to the Jeans mass criterion it can spawn a 
star particle with a mass one eighth of the mass of a gas particle (thus a given gas particle can 
produce at most eight star particles). The process is governed by either drawing from a Poisson 
distribution such that the local rate of star formation agrees with Eq. 19 in a stochastic sense, or 
by imposing a fixed delay time. 

Our basic star formation recipe does not depend explicitly on the local H2 fraction, so while 
the gravitational collapse that induces star formation may also precipitate H2 formation, it is 
nevertheless independent of whether the gas is atomic or molecular. However, our model now 
allows for a direct link between star formation and the presence of molecular gas. We can do this 
by setting a local threshold for fm above which star formation can proceed. As we will discuss in 
section 4.2.2 this may be more realistic since in such a scenario the values of the delay factor are 
no longer assumed a priori but are instead implemented in a more physical fashion by the "delay" 
associated with the final chemical/ thermodynamic evolution of a gas cloud before star formation. 

In order to determine the local FUV field used for calculating the photoelectric heating and 
the H2 destruction rates, we determine the time-dependent FUV luminosity of the stellar particles. 
We do so by following their age and, since star particles represent stellar associations rather than 
single stars, by using Bruzual & Chariot (1993 and updated) population synthesis models. We 
assume a Salpeter initial mass function (IMF) with cutoffs at O.IMq and IOOMq. In the present 
work we do not account for dust extinction of UV light, except that of the young stars shrouded in 
their natal cloud: for a young stellar cluster we decrease the amount of UV extinction from 75% 
to 0% in 4 Myr (Parravano et al. 2003). 

Feedback from stellar winds and supernovae is essential for regulating the physical conditions 
of the ISM. While the mechanical energy output of stars is reasonably well known, it has proven 
to be difficult to include it completely self-consistently in galaxy-sized simulations of the ISM. The 
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reason for this is that the effective energy of feedback depends sensitively on the energy radiated 
away in thin shells around the bubbles created. This will mean that the effect of feedback cannot 
be tracked in a straightforward manner unless prohibitively high resolution is used. In SPH codes 
there have been conventionally two ways to account for this: by changing the thermal energy 
input and by acting on particle velocities. Both are unsatisfactory, as the thermal method suffers 
from overcooling (Katz 1992) and the kinetic method seems to be too efficient in stirring the ISM 
(Navarro & White 1993). Here we use a new method based on the creation at the site of young 
stellar clusters of a pressure particle that acts as a normal SPH particle in the limit of the mass 
of the particle m ^ 0, for constant energy. For the energy injection rate (which assumes that 
the energy injection takes place continuously, without distinguishing between stellar winds and 
supernovae) we take 

E = esnHsnEsn/At, (20) 

with Egn = 10^^ erg the energy liberated per supernova, an efficiency parameter Cgn = 0.1, the 
number of supernova per mass of new stars formed Ugn = 0.009 per Mq (appropriate for the 
IMF adopted here) and At = 3 x 10^ yr. The efficiency egn thus assumes that 90% of the ini- 
tial supernova energy is radiated away in structures not resolved by our simulation. This value 
has been found in detailed simulations of the effects of supernovae and stellar winds on the ISM 
(Silich et al. 1996; Thornton et al. 1998), and is also used in other simulations of galaxy evolution 
(Semelin &; Combes 2002; Springel & Hernquist 2002; Buonomo et al. 2000). Moreover in a study 
by Pelupessy et al. (2004) the adopted feedback strength has been shown to be consistent with the 
observed scatter in the star formation properties of isolated dwarf galaxies. Details of the feedback 
model, including implementation features and observational constraints can be found in Pelupessy 
et al. (2004) and Pelupessy (2005). The most important effect of SN feedback for our model of 
H2 formation is, apart from the regulating effect of feedback on star formation, the introduction of 
extra dynamical pressure (Eq. 11) through the increases in local velocity dispersion. 



3.3. H2 formation 

For our H2 formation model we follow the same philosophy as for the star formation model: 
unresolved structure is assumed to be present at the SPH particle positions. In this case the 
underlying structures are assumed to conform to Eq. 10 expressing the mean extinction and 
derived from the observed density/size scaling relation (Eq. 8). However the latter is not expected 
to be valid throughout our simulation domain. Regions where the density/size relation becomes 
inapplicable are those with low density and pressure where it predicts very large cloud sizes and the 
resulting photo-destruction of H2 proceeds very slowly. Such regions contain WNM HI gas where 
the pressure is mostly thermal and then Eq. 8 (along with Pe = UekeTk) yields Rpc oc T\^Pe^^^, 
which for typical WNM conditions corresponds to kpc-size "clouds" . To circumvent this problem 
we modify Eq. 8 for Pe < Ptrans as: 



where we take Ptrans ~ 1000 K cm~^. This is just a convenient patch, which however in low 
pressure environments does scale the "cloud" mean density as oc Pc, found for diffuse clouds 
(Elmcgrecn 1993). This minor modification now allows the application of our subgrid model in 
all ISM conditions present in the simulations. 

For the macroscopic pressure Pe in Eq. 11 we need the local velocity dispersion cTyei- For this 
we take the formal SPH estimate 



with Vi and mj the particle velocities and masses, (v)j the local bulk velocity. A number of subtle 
issues are connected with the choice of the dispersion. Equation 22 describes the inter-particle 
velocity dispersion, while the dispersion that enters Eq. 11 is really an intercloud velocity dispersion. 
Even if we assume these to be equivalent, the problem remains that Eq. 22 expresses a velocity 
dispersion on scales of the local SPH smoothing length h^. We can try to account for this by scaling 
(jj, using for example the relation for Kolmogorov turbulence. 



but we found this to have rather little influence on the resulting H2 formation. 

Depending on the cloud model, we use Eqs 12 or B5 when T^ < 1000 K to track the evolution 

(tr) 

of Av , and thus fm, during a simulation timestep dt. For T^ > 1000 K and Rf = Eq. A3 
describes the photodestruction of clouds, while for T^ > 3000 K, Eq. 17 is used to approximate 
the coUisional destruction process of the remnant molecular gas. The density that enters those 
equations is assumed to be the mean density (n) given by the SPH density at the particle position, 
and the temperature the particle temperature (both taken constant for the timestep) . The radiation 
field is calculated from the distribution of stars, where the assumption is that extinction from dust, 
or from molecular clouds is not important (apart from extinction from the natal cloud). After we 

(tr) 

have evolved Av for the timestep dt the resulting fm (Eqs 7, B3) is calculated and assigned to the 
SPH particle for the duration of the next timestep (where it can be used to calculate e.g. cooling). 
At the next timestep the last fm value is retained and used to calculate the initial Af) , given the 
new average cloud extinction (Ay), which may have changed from the previous time-step, if for 
example the pressure Pe has changed. 

The choice to keep fm rather than A^^^ constant during any variations of (Ay) is dictated 
by the decision to assign all variations of fm to the FUV-regulated HI/H2 gas phase interplay, as 




(22) 




(23) 
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tracked by Eqs 12 and B5 rather than to any instantaneous effects. Indeed, if instead of 

fm was to be kept constant during variations of (A^) over a given time-step, it would correspond 
to an fjn that instantaneously follows the variations of (Av). This may not be entirely without 
merit since whatever processes are responsible for setting up cloud structures that obey power 
laws like that expressed by Eq. 8 may be also responsible for a fast concurrent H2 formation 
(Chieze & Des Forets 1987; De Boisanger k, Chieze 1991). Nevertheless at this stage we chose not 
to consider this possibility since the cloud power laws are assumed "frozen" during our simulations 
and allowing instantaneous effects on fm (via variations of (Av)) would run counter to the spirit 
of our effort to model the time evolution of the HI, H2 gas phases without such influences. In the 
future we intend to explore this possibility to see whether it is important when compared to the 
FUV-regulated HI/H2 phase interplay explored here. 



3.3.1. The applicability of the H2 formation model at high resolutions 

In order to track dense gas our simulations employ SPH particle masses of mspH = 500 M© 
(see section 4.2.1) where it can be argued that the applicability of the scaling laws used to deduce 
fm is doubtful because such small masses are well below those of even the smallest GMCs and thus 
represent cloud fragments instead. However, individual SPH particles do not represent the smallest 
resolved objects of our simulations, rather this is set by the total number of particles used to derive 
the formal SPH estimate of the mean density. The latter is of the order Mj-ef ~ 10^ Mq (see Section 
4.2), comparable to small molecular clouds. Furthermore, it must be noted that the applicability 
of the n-R scaling law for their Mj-ef sub-units is based on the virial theorem and the universality of 
the linewidth-size relation. The latter has now been confirmed to hold for scales and environments 
well inside typical molecular clouds and with the same normalizing constant as for entire clouds, a 
fact attributed to the universality of large-scale driving mechanisms for turbulence (Heyer &: Brunt 
2004). Regarding the virial theorem, it is a simple corollary that if GMCs are virialized objects so 
will be any of their sub-units albeit at different boundary pressures, thus Eqs 8, 10 are expected to 
hold even for M = Mref cloud masses making up larger GMC-type ensembles (~ (10^ — 1O^)M0). 
The latter of course assumes that the correct boundary pressure Pe is used in Eqs 8, 10, which as 
we already discussed may be poorly represented by the formal SPH estimate of the macroscopic 
velocity dispersion. Thus our model is applicable for scales where the size-linewidth relation and 
the assumption of virial equilibrium are valid, which will certainly not be the case in dense star 
forming cores, but does seem to be true at the larger scales as probed by galaxy scale simulations. 



4. An application to dwarf galaxies 

Our first application will be a model of a dwarf irregular galaxy. There are a number of practical 
reasons to choose this type of a system as a test model, as well as some interesting issues that are 
particular for dwarf galaxies that can be explored using our model (e.g. H2 vs HI gas supply) . The 



-16- 



small size of these systems allow relatively high resolution with modest computational effort. Indeed 
the choice of dwarf galaxies as modeling templates allows a given numerical simulation to probe 
small physical scales and high gas densities, the latter being of crucial importance if the HI II2 
phase transition is to be tracked successfully. Insight gained from simulating these systems, whose 
properties are constrained by a wealth of observational data (Van Zee 2001; Barone et al. 2000; 
de Paz, Madore & Pevunova 2003). Furthermore it holds the promise of yielding good constraints 
about physical processes that are expected to be universal in galaxies (see e.g. similar work done 
for SN feedback strength; Pelupessy, van der Werf, &; Icke 2004). 

Aside from being excellent testbeds for quantifying phenomena expected to be common across 
the Hubble sequence, dwarf galaxies are important systems on their own right because of their 
possibly very important role in current galaxy formation theories as the building "blocks" of larger 
systems (Kauffmann et al. 1993), and as major "polluters" of the intergalactic medium with metals 
(Ferrara & Tolstoy 2000). 



4.1. Simulation setup 

We construct a simple model of a dwarf irregular with gas mass of Mgas = 10^ Mq and a stellar 
mass of Mgtar = 10^ Mq. Specifically, the gas disk has radial surface density profile 

s = Sg/(i + R/Rg)> (24) 

with central density Sg = 10 Mq/pc^ and radial scale Rg = 0.33 kpc, truncated at 4 kpc. An 
exponential stellar disk, 

Pdisk(R, = ^ exp(-R/Rd)sech2 (z/h,) (25) 

with central surface density Sq = 300 Mq/pc^, Rd = 0.5 kpc and vertical scale height h^ = 0.2 kpc, 
is constructed as in Kuijken & Dubinski (1995). The gaseous and stellar components are repre- 
sented by 2 X 10^ particles each. For the gravitational softening length of the stellar particles a 
value of 20 pc is adopted, while for the gas particles it is taken to be equal to the SPH smoothing 
length. The ages of the initial population of stars are distributed according to a constant star 
formation rate (SFR) and an age of 13 Myr. 

The rotation curves of dwarf galaxies are generally best fit using dark halos with a flat central 
core (Flores &; Primack 1994, Burkert 1995). Therefore we take a halo profile 

Phalo(r)=po 7|7//? (26) 
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with core radius 7 = 2 kpc, cutoff radius rc = 20 kpc and central density po = 2 x 10^ M0/kpc^, 
for a total mass of Mhaio = 15 x 10^ and a peak rotation velocity of about 50 km/s. Note 
that the profile adopted is very similar to the Burkert (1995) profile. They differ mainly in their 
asymptotic behaviour for r ^ 00, thus the Burkert profile will only deviate significantly well outside 
the region of interest for our simulation (within 20 kpc the difference in rotation velocity of the 
adopted density profile and a Burkert profile with the same central density is < 6%). We represent 
the dark-matter halo by a static potential. 

4.2. Practical considerations 

The mass resolution of the simulation limits the density that can be probed. Wc can derive 
some minimum requirements in order for our model to be able to follow the HI H2 transition. 
In addition, it may be necessary to consider the effect of II2 cooling. Finally, unlike previous 
numerical work reported in the literature, the range of values for the delay factor fgf , which regulates 
star formation, can now be explored under constraints imposed by the observationally-motivated 
demand that H2 forms ahead of stars. 

4-2.1. The necessary resolution: the choice of Mref 

Simulations of self-gravitating fluids done with insufficient numerical resolution can suffer from 

artifacts: artificial clumping or inhibition of gravitational collapse may occur. For SPH simulations 
including self-gravity these effects can arise if the local Jeans mass is not resolved (Bate &; Burkert 
1997; Whitworth 1998). Hence, this requires for our simulation 

Mj > NmsPH, (27) 

with N the number of SPH neighbours and mspH is the mass of an SPH particle. In our simulations 
this requirement is met by choosing appropriate star formation parameters so that violation of Eq. 
27 is precluded by star formation. A gas particle will "spawn" star particles, and thus be subjected 
to heating that raises the local Jeans mass, whenever Mj < Mj-cf- Taking Mrcf = NmspH in effect 
makes the choice of mass (and density) resolution equivalent to the choice of Mref- The density 
resolution of our simulations must be high enough so that H2 formation can compete with FUV 
dissociation (i.e. T^is < !> Eq. 13). In Figure 1 it is apparent that densities of ~ (50 — 100) cm~"^ are 
sufficient to follow the HI H2 transition, while somewhat higher densities may be necessary for low 
metaUicity gas. For such densities and at typical Tk(CNM) ~ 100 K, Eq. 18 yields Mj ?s 10^ Mq. 
Hence running the simulation with particle masses of mspH ~ 500 Mq allows a good resolution of 
the aforementioned mass scale. Apart from properly tracking the HI/H2 phase transition, setting 
a small Mref ~ 10^ allows GMC-type, star-forming associations of gas particles to naturally 
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emerge in the simulations. Indeed, GMC-class gas masses of Mgmc ~ (10^ - 10^) M© can then 
emerge as asscmbhes of smaller star- forming cloud units. Our star formation criterion Mj < Mj-cf, 
where Mj-cf < MqmC) ensures that such self-gravitating SPH particle formations (which previous 
simulations show emerging, e.g. Gerritsen 1997) always contain gravitationally unstable regions 
that form stars, exactly as in real GMCs. In that respect Mj-ef can be viewed as a mass "radius" 
attached to individual SPH particles which defines the Jeans-mass instability, and which is itself 
smaller than the mass of typical GMC aggregates, the preferred sites of star formation in galaxies. 
If one were to choose Mj-cf ^ Mgmc stars could then form from "cloud" mass scales that are 
never seen star forming in nature. As a final note on the choice of Mj-cf it worths mentioning 
that observational studies have shown that O, B stars do not form in GMCs with masses below 
~ 10^ Mq (Elmegreen 1990). Thus choosing Mref ~ 10^ M© makes our simulations capable of 
following the evolution of the smallest GMCs seen forming stars in galaxies. 



In our standard star formation recipe the delay for star formation to occur is set by the free 
parameter fgf; our default choice is a fixed value of fgf = 10. In the past typical values of fgf ~ 5 — 10 
were chosen so that the expected star formation rate from the GMCs present e.g. in the Milky 
Way would be similar to the one observed. Dedicated high-resolution simulations of individual 
molecular clouds may shed some light on this issue (e.g. MacLow et al 1998), but for the time 
being the choice of fgf in galaxy-size simulations like ours remains rather apriori. 

Nevertheless the fact that everywhere in the local Universe stars seem to form out of molecular 
rather than atomic gas allows us to deduce a lower limit on the value of fgf. Indeed, for H2 formation 
to remain always ahead of star formation it must be rf(H2) < Tgf, which from Eqs 2 and 19 yields. 



Since CNM HI is the most likely precursor phase of the H2 gas, for typical conditions of Tk ~ 
(100 — 200) K and (ni) ~ (10 — 50)cm~^, the latter equation yields fgf > 3 — 10. These values are 
roughly similar to those deduced using the constraints on the observed star formation rate in the 
Galaxy. Given the uncertainties inherent in e.g. the SH(Tk) function, such a rough agreement is 
noteworthy and here it is worth mentioning that rf(H2) is also a good approximation of chemical 
equilibrium timescales in FUV- illuminated clouds (Hollenbach & Tielens 1999). Thus the rough 
agreement of the constraint set by Eq. 28 with what are considered reasonable fgf values may 
signify an important role for cloud chemistry in setting the average star formation rate. This is not 
far-fetched since molecule formation enables very different gas cooling functions and thus drastically 
alters its thermodynamic state allowing to eventually cool and fragment further (well beyond the 
resolution limit for galaxy-sized simulations) (e.g. Chieze &: des Forets 1987). For these reasons we 
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will also discuss a model for the star formation which depends on fm directly. 

4.2.3. Cooling by H2 

A detailed calculation of H2 cooling was done by Le Bourlot et al. (1999). At Tk ~ (1000 — 
8000)K H2 is a more efficient coolant than the other major coolants (C and Fe). Hence even a 
low (~ 10^^) remnant abundance of molecular gas in the diffuse WNM can have an impact on 
the cooling at these temperatures. Indeed some people have included H2 cooling in their cooling 
curves, but without calculating the abundance of molecules (Carraro et al. 1998). Admittedly our 
accuracy in following the amount of warm molecular gas is limited, but to get some idea of the 
possible effect of such a gas phase it is nevertheless interesting to do some calculations including 
H2 cooling. For this we use the Le Bourlot cooling curve, where we use a limited set of data 
(low density limit, constant ortho-to-para ratio, only H2-H collision excitation), appropriate for 
our purposes. Note that some observations find abundances of H2 gas the order of 10~^ — 10~^ in 
diffuse HI gas, even in regions of the ISM hostile to H2 formation, like galactic high velocity clouds 
(Richter et al. 2001). 

5. Dwarf galaxy model: results 

In a model that incorporates H2 formation the two obviously important parameters that should 
be explored are: the formation rate /x, and the mctallicity Z. The formation rate is uncertain by a 
factor of 5-10, and thus introduces some uncertainty in the theoretically calculated H2 content, e.g. 
for high formation rates significant amounts of warm and diffuse H2 may be present in the outer 
parts of spiral galaxies (Papadopoulos, Thi, &: Viti 2002). Metals play an important role: the H2 
formation rate is directly proportional to the dust surface available, and we take this surface to be 
proportional to Z. Dust also plays an important role shielding molecules from UV radiation. Obser- 
vations of low metallicity systems, like dwarf irregulars, using CO as a tracer give low molecular gas 
contents (Barone et al. 2000). However, the interpretation of these observations is complicated by 
the fact that the conversion factor from CO flux to H2 mass under these conditions depends heavily 
on metallicity and ambient FUV field: for low metallicity, FUV-intense environments CO dissociates 
while the largely self-shielding H2 is not affected by much (Maloney & Black 1988; Israel 1997). 

In this work we explore a 2x2 grid of models using a formation rate of = 3.5 and ji = 17.5, 
and metallicities of either Z = Z0/5 or Z = Z©. Note that for our model dwarf galaxy the high 
metallicity is somewhat unrealistic, because normally dwarf galaxies have lower metal content. Note 
also that for simulations without H2 cooling our subgrid model follows the H2 content as a passive 
tracer of the history of the local conditions, thus the simulations for high and low formation rate 
will be identical except for the amount of H2 . 
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5.1. Star formation properties of the simulations 

We start the simulation with an isothermal gas disk at ~ 8000 K, and a molecular fraction 
fm = 0. After the start of the simulation the gas cools and collapses along the z-direction, forming 
regions of cold and dense gas in the midplane (the minimum size of structures resolved in the 
simulation is about 30 pc). There, H2 starts to form and star formation starts taking place as well. 
The fact that we start from an unrealistically smooth initial condition results in some transient 
effects, but after approximately 200 Myr the galaxies settle in a stationary state, with a constant 
star formation (see Figure 4). For this model galaxy, the mean star formation rate is about 
SFR = 0.002 Mg/yr for low Z and SFR = 0.003 Mq/jt for high Z, with modest variations (the gas 
depletion timescale is much longer than the duration of our simulation, in the order of 5 x 10^° yr). 
The star formation values wc get represent fairly typical values for a low surface brightness dwarf 
(Van Zee 2001). Also plotted in Figure 4 is the average star formation density as a function of 
radius, note the exponential distribution of star formation and the sudden truncation. During the 
simulation a small fraction (1 x 10^^) of the gas is expelled out of the disk, but most of this falls 
back. A hot galactic wind is not formed, the low mass loss rate associated with the expected wind 
for this model cannot be resolved by our current mass resolution. 

5.2. Spatial and temporal distribution of H2 gas 

An equilibrium fraction of H2 is reached on timescales of 50 — 100 Myr (see Figure 5), 
in agreement with our earlier approximate estimates from Eq. 2. Furthermore, as expected, the 
low-/Lt/low-Z simulation yields the lowest equilibrium molecular fraction of fm 0.001 and a higher 
formation rate n = 17.5 will boost that to fm ~ 0.03. For high metallicities, molecular fractions of 
fm ~ 0.17 {iJ, = 3.5) and fm ~ 0.4 {jj, = 17.5) form. Also drawn in Figure 5 is the time dependence of 
the molecular fraction for a low-/i/low-Z simulation with logotropic clouds. The difference between 
constant density and logotropic clouds is not very large, certainly smaller than the differences due 
to the uncertainties of ^ (this is also evident in the equilibrium values of fm shown in the panels in 
Figure 1, where the main difference between constant density and logotropic clouds was that for 
logotropic clouds appreciable fractions fm appear also at lower densities). 

In Figures 6 and 7 we show the (mean) molecular fraction as a function of radius and height 
above the disk-plane for the four simulations. We see that, as expected, H2 forms mostly in the 
central regions and in the midplane where the pressure and gas density builds up. Note also that 
for high metallicity and high formation rate the molecular fraction seems to have a tendency to 
"saturate" in the central regions at a fixed value, in this case fm = 0.6. Looking in more detail 
at the spatial distribution of HI and II2 shown in Figure 8 we find II2 mainly concentrated in the 
dense clumps and filaments of the HI distribution. The structure of the ISM can be quantitatively 
compared with observations by analyzing the power spectrum of the HI maps such as those in 
Figure 8. This is done in Figure 9, where we see that it agrees quite well with power spectra 
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made for the most detailed studied dwarf galaxy, the Large Magellanic Cloud (LMC), and is stable 
over the course of the simulation and as a function of resolution. If we now compare the low and 
high fi simulations, we see similar H2 distributions: the uncertainty in the formation parameter 
is mainly an uncertainty in the amount of H2 formed, not in the locations of H2 formation. The 
relative distribution of H2 is well represented by our simulation, but the exact amount of H2 is 
more difficult to determine, as this depends on the formation rate parameter. If we look at the high 
metallicity simulations, wc can sec that, apart from the clumpy distribution there is also a sizable 
mass of "diffuse" H2, but still confined mainly in gas filaments. 

Finally we look at the relation between star formation and H2 distribution. In Figure 10 
we show a simulated "Ha" map (actually a map of the stellar luminosity in ionizing photons), 
overplotted with contours of H2 gas, chosen so that they highlight the relation of star formation 
to the densest parts of the H2 distribution. It can be seen that star formation is almost always 
associated with some nearby H2 gas complex while some of them, being without newly formed stars 
at that particular instant, do not show any Ha emission. 

5.3. The influence of the collapse delay factor fgf 

Apart from the parameters directly related to H2 formation, fm will also depend on the delay 
factor fsf. As discussed in section 4.2.2, this parameter accounts for the fact that a region that is 
Jeans- unstable does not collapse to form stars at the free-fall time tq. The total amount of H2 
is sensitive to the time HI gas stays in a cold and dense phase conducive to H2 formation, while 
the emergence of new stars from such a phase helps dissociate H2 by dramatically increasing the 
ambient FUV radiation. Indeed, comparing models run with different values for fgf reveals the 
total amount of molecular gas to be strongly dependent on this parameter (Figure 11). In short, 
the formation of stars and H2 arc in competition, and small enough delay factors can quench H2 
formation. Yet, a value of fgf > 10, as constrained by the observed GMC masses and the star 
formation rate in the Galaxy, amply satisfy Eq. 28 and H2 formation preceeds star formation. 

5.4. Star formation with an H2 gas fraction threshold 

In section 4.2.2 we have argued for a connection between star formation and chemical tinicscalcs, 
the latter well approximated by the H2 formation timescale. A simple and direct implementation of 
this idea is to replace the preset delay time Tgf by a threshold molecular fraction fm^gf for the onset 
of star formation. In this case the delay factor fgf is no longer needed; the "delay" from a pure 
free-fall timescale is now the outcome of the interplay between the competing physical processes 
that lie behind the establishment of fm- 

We rerun the Z = Zqj/x = 3.5 simulation with this new star formation recipe and a threshold 
value of fm^sf = 1/8. In figure 11 we have plotted fm for this case. As can be seen, the equilibrium 
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molecular fraction goes down significantly, to fm ~ 0.01, which is easily understood, because such a 

low threshold value will almost immediately destroy H2 once it forms. For higher threshold values 
fm is higher, fm ~ 0.05 and 0.1 for fm,sf = 3/8 and 5/8 respectively (this is still less than for 
the normal recipe, because for this case about 40% of H2 is in regions with fm > 5/8). The star 
formation rate for these models is about SFR = 0.004 Mq/jt. 

In figure 12 we plot the resulting (cumulative) distribution of delay times Tgf, normalized 
on the local free-fall time tg. A higher threshold results in longer average delay times, from 
{Tsf) ~ 1.5tff for fin,sf = 1/8 to (rgf) ~ lOtg for fm.sf = 5/8. However, the fraction with Tgf < tff 
stays roughly constant. The reason for this is that once star formation starts in a region, it takes 
some time for feedback to kick in. In the mean time molecule formation continues, and the molecular 
regulated star formation allows for extra star formation events if enough H2 forms (in the previous 
implementation of star formation there was no such possibility). In other words, the local star 
formation efficiency is partly determined by the H2 formation. 

A high value for fm,sf seems to be favoured, because star formation regions are observed to be 
mainly molecular, and in this case both the H2 fractions and local star formation rates are consistent 
with what is known from observations. Hence this model for star formation seems promising for 
future application. 

5.5. Simulations with H2 cooling 

The inclusion of H2 cooling affects the simulation in a number of ways. As an extra coolant it 
will increase the amount of gas in a cold state, and in turn this will increase the formation of H2. 
However, star formation will also increase, and thus also the UV field. In Figure 13 the effect of 
cooling on fm is illustrated for the ji = 3.5 simulations. As can be seen, cooling has some influence 
for the high Z simulation but for low Z the molecular fractions are too low to substantially alter the 
results. Note that while H2 is a strong coolant at high temperatures (a few thousand K), where it 
dominates the cooling for fm > 0.001, the effect of the extra cooling is limited because of the small 
amount of gas at these temperatures. 

6. Discussion 

Detection and mapping of the molecular component of dwarf galaxies is usually done by map- 
ping the ^^CO(J = 1 — 0) line (Barone et al. 2000; Mizuno ct al. 2001) or by UV absorption studies 
(Tumlinson et al. 2002). The derived H2 fractions are ~ (1 — 10)%, but for low metallicity systems 
CO is often not detected corresponding to upper limits for the molecular fraction of ~ 1%. Recently 
a large observational effort has detected ^^CO J=l-0 in several more dwarf galaxies increasing the 
number of such systems with detected CO emission by ~ 50% (Leroy et al. 2005). 
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At this stage a direct comparison of our simulations with available CO imaging data must be 
made cautiously. This is because our simulations follow the formation of H2 not CO. The latter is 
the most abundant molecule after H2 itself and serves as the prime tracer of its mass (e.g. Dickman, 
Snell, & Schloerb 1986; Young & Scoville 1991), yet it is still four orders of magnitude less abundant 
and thus, unlike H2, it cannot self-shield. This is the main reason why in metal-poor and FUV- 
intense environments like those prevailing in dwarf galaxies CO (but not H2) will be dissociated 
and its emission can then give a rather misleading picture of the H2 mass distribution (e.g. Madden 
et al. 1997). Even if a simplified version of the chemical network giving rise to CO (e.g. Nelson 
& Langer 1997) were to be included in the simulations, densities n > 500 cm^"^ and temperatures 
Tk < 50 K would have to be tracked in order to successfully model CO formation, and these are 
currently inaccessible by galaxy-size simulations like ours. 

In principle one could use empirical relations of the X = M(H2)/Lco(i_o) = F(Go,Z,n, Tk) 
factor available from the literature (e.g. Israel 1997; Bryant &; Scoville 1996), and then convert 
the H2 maps from our simulations to CO brightness maps. This we intent to do in future work 
that will also include spiral disks. Nonetheless, the molecular fractions we derive here seem to 
be reasonable, although for the low-/i/low-Z simulations the derived fm is on the low side. The 
molecular fraction may increase somewhat for a simulation run at higher resolution, because in 
metal-poor environments the HI H2 transition is taking place above n = 100 cm~^, this can 
increase fm by a factor 2-5. For the galaxy models run with high metallicity and/or formation rate 
we find quite substantial fractions of H2, of fm = 0.17 — 0.4. 

We also seem to find quite high values, fnj ~ 0.001 — 0.01, of the diffuse warm neutral medium 
to be molecular. To illustrate this we have plotted in Figure 14 the molecular fraction as a function 
of total column density for a face-on projection. There is a clear column density threshold for the 

presence of H2 for the low metallicity simulation. For high metallicity even low column density 
regions (which arc also of low density and high temperature) show fractions of fm ~ 0.001 — 0.01. 
This state of the gas (low density, high temperature) is not typically associated with the presence 
of such quantities of H2. There is no - or very little - formation of molecular gas going on in 
this gas phase and any H2 present should dissociate. In fact, destruction is taking place in the 
simulation, but its timescales are quite long. For the thermal dissociation this is due to the low 
density of the gas: for temperatures of T = 5000 K and density of uh = 0.3 cm~'^ the collisional 
destruction timescale (for neutral gas) is TcoI ~ 10^ yr. Furthermore, for our model galaxy the 
radiation fields in quiescent regions are quite low. Go ^ 0.1, due to the low star formation, so 
the radiative dissociation timescales are also in the order of Trad = 1.3/Go x 10^ yr 10^ yr. 
Hence the H2 survival in lukewarm gas may not be unrealistic. Note however that the simulation 
may not represent the dispersion of H2 from dense regions realistically: molecular gas is formed in 
dense, star forming regions, where radiation fields are higher (for our model up to Go = 100. As 
the star formation heats the surrounding gas it will destroy H2, the amount of molecular gas that 
"escapes" destruction is then sensitive to e.g. the time it is exposed to the high radiation field or 
the combination of high density/ high temperature. Also, destruction by supernova shock is not 
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represented very well, because the interaction of the supernova shocks with the (sub resolution) 
molecular clouds is absent. On the other hand shocks have been shown recently to also induce fast 
H2 formation as they propagate through diffuse HI gas (Bergin ct al. 2004) and the effect of such 
a mechanism in the overall HI/H2 balance in our simulations is unknown. 

Nevertheless, our simulations do indicate that in so far as the molecular gas is not destroyed 
at the sites of its original assembly, it can survive the CNM-to-WNM transition. In that case 
the lukewarm and WNM gas phases may still contain significant amounts of H2. In subsequent 
work we intend to explore the issue of H2 survivability under such conditions more thoroughly by 
introducing a better model of its coUisional and FUV-induced destruction in such gas phases while 
also including shock-induced H2 formation and destruction mechanisms. Finally, any significant 
underlying cloud substructure in CNM HI clouds is likely to enhance the amount of H2 gas present 
in that phase. 

Future applications of H2-tracking numerical models of galaxies can be divided in two broad 
and complementary categories, namely: a) examine the effects of the various model parameters 
on the H2 distribution within a given galaxy by e.g. using different H2 formation functions (see 
Cazaux & Tielens 2004), and/or using a range of fgf, esF and Mj-ef values, b) explore the evolution 
of the H2/HI distribution across the Hubble sequence using a single "mean" H2-tracking model. In 
the latter category starbursts/mergers, with their extreme and fast-evolving pressure and radiation 
environments (and hence probably strongly varying H2/HI fractions) provide particularly attractive 
templates for modeling. 

7. Conclusions 

We presented a method to calculate the local H2 content in simulations of the ISM based 
on a subgrid model for the formation and destruction of molecular gas. The model tracks the 
formation of H2 on dust grains and its destruction by UV irradiation in the CNM phase and 
its collisional destruction in the WNM phase including the effects of shielding by dust and H2 
self-shielding. The solution to the H2 formation/destruction problem is simplified greatly by the 
assumption that H2 formation takes place in structures that conform to the observed density /size 
scaling relations, present on scales not resolved by the hydrodynamic code of the simulation. We 
show steady state solutions of the model for the molecular fraction fm as a function of the average 
gas density, temperature, velocity dispersion, radiation field, and metallicity. The effects of the 
cloud substructure are explored by taking 1/r density profiles for the model clouds, this seem to 
have only a modest effect on the derived H2 fractions, but it does allow more diffuse gas to be 
molecular. 

A simple time-dependent formulation of the same model is then coupled to hydrodynamic 
simulations and used to calculate the local evolution of the molecular gas fraction according to 
the local macroscopic quantities. We have then incorporated our H2-tracking method in an N- 
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body/SPH code for the simulation of galaxy-sized objects. Our first application is to model a typical 

low surface brightness dwarf galaxy where the demands set by the density resolution necessary to 
properly track the HI — > H2 transition are easily met. Our results can be summarized as follows: 

• We found that our model reproduces reasonable molecular fractions ranging from fm = 0.001 — 
0.03 for low metallicity to fm = 0.17 — 0.4 for solar metallicity systems, saturating to fm = 0.6 
for high values of the H2 formation rate and in the central regions of our model galaxies. 

• The biggest uncertainty is the value of the adopted formation rate parameter n, however the 
general features of the spatial H2 distribution do not depend significantly on /i. 

• The molecular fraction shows a strong dependence on the metallicity Z, as expected for H2 
forming on dust grains. Our model for a high metallicity dwarf shows surprising amounts of 
"diffuse" H2, possibly due to the low radiation fields outside star forming regions. This may 
be analogous to what happens in the outer regions of spiral galaxies where the radiation field 
drops, and a diffuse H2 gas phase may be present there. However a better model for the H2 
formation and destruction in the WNM phase is needed before the aforementioned results 
can be considered secure. 

• The molecular gas content is sensitive to the details of our star formation recipe, especially 
on the value of the delay factor fgf adopted. This poorly constrained factor, generally used 
to parameterize a slower-than-free fall star formation timescale, needs to be sufficiently large, 
fsf ^ 10 to allow for H2 formation before young stars start to dissociate H2. On the other 
hand a star formation criterion based on the formation of H2 can dispense with the fgf factor 
completely, giving local star formation rates well within constraints set by observations. This 
suggests that the star formation is governed by cloud chemistry processes, and indicates the 
importance of considering molecular gas in simulations. 

• H2 cooling has only a minor impact for our model. The largest effect of H2 cooling is seen in 
high metalicity (and thus H2 content) environments. 

Summarizing, we developed an algorithm to calculate the evolution of the molecular gas phase 
during the evolution of the gaseous and stellar content of galaxies, and obtained a first application 
to dwarf galaxies. Strong points are that it is simple, physically motivated and not computationally 
expensive. It is well suited to galaxy-sized simulations where a complete treatment of the physics 
involved in the formation of the H2 clouds is inaccessible, while our approach can be seen as 
complementary of (much higher resolution) simulations modeling individual molecular clouds. The 
set of astrophysical questions that can be addressed using such H2-tracking models is large. These 
include examining the influence of "micro" -physics parameters/functions like the H2 formation rate 
on the general molecular gas distribution within a given galaxy, as well as modeling the gaseous 
ISM evolution for galaxies across the Hubble sequence. The latter can now be done by including 
all gas phases along with the stars, with the molecular phase expected to be the most intimately 
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involved with the star formation. Finally, mergers/star bursts, with their fast-evolving pressure and 
radiation environments are ideal templates for future applications of the model presented here. 
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A. Analytical approximations for A^*'^) 

It is interesting to examine approximate analytical solutions of Eq. 12, if only for the reason of 
checking the more general solution in limiting cases. This can be done for two particular domains, 
namely that of rapidly increasing and rapidly decreasing Apy^y(t). For TfdAp^y/dt ^ we can 
approximate Eq. 12 as 

— ^^-f^FUV ^ ^FUV ~ ^FUV*.^/^ • \-^^) 

In the case of TfdAp^y/dt S> 1, and thus a decreasing molecular mass fraction, we approximate 
Eq. 12 as 

= ^e-<V = 2Goko$e-<V, (A2) 

dt Tf o o , y J 

hence 

AtuV(t) = In (e4uV(o) + 2Goko$t) , (A3) 

which, for Rf = 0, provides a convenient solution describing the photo-destruction of molecular 
clouds. 



B. logotropic clouds 

Here we give the corresponding equations of section 2 for clouds with a density profile n(r) = 
ne(r/R)~^ ("logotropic clouds"). For a radially incident interstellar radiation field, such a density 
profile yields 

Ntr(HI) = ^lnfl + ^-i^$y (Bl) 

where i^o = UeRcj (1 + UeRcr)^^. For R ^ oo we obtain i/q ^ 1 and the transition column density 
of Eq. Bl reduces that of Eq. 4 (for n(r) Ue), as expected, because for large clouds the density 
will not change much over Ntr(HI). The visual extinction corresponding to this Ntr(HI) is given by 

Af) = 1.086i/oCpuvln 



/xSH(Tk) V ZTk J 



135 cm ^ 



(B2) 
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and the H2 fraction fm, 



M(H2) 



exp 



(Av) 



(B3) 



The modification of Eq. 12 for the time dependent transition column in the case of a logotropic 
density profile (Eq. 22 of Goldshmidt & Sternberg 1995) turns out to be particularly simple, namely 



ldN(HI,t) Goko$ _^ 



N(HI,t) 



dt 



Rf / n(r)n(HI)dr. 



(B4) 



Here rni marks the depth of the HI layer. Using the relation n(r) = Uee^^''-*/^" (No = UeR) valid 
for a logotropic density profile, Eq. B4 eventually yields 



Tf- 



dCTNu(Ill,t) 
"dt 



rdise-"N-("i'*) - ctNo (e'^N,.(Hi,t)/<.No _ ^ 



(B5) 



where Tf and r^is are calculated for n = ne. It can be easily seen that for No = n^R — >^ 00 B5 
reverts back to Eq. 12 with a density n = ne (i.e. the case of a cloud with so large a radius that the 
HI/H2 equilibrium is well approximated by a plane-parallel geometry at a uniform density n = ne). 

ftr) 

For dA^{/y(t)/dt » the solution of Eq. B5 remains the same as in Eq. 12 but with 
ne = 2/3(n) replacing n. Hence, in the case of cloud destruction and Rf = 0, Eq. A3 expresses the 
time dependence of the transition layer also in logotropic clouds. For dApuy(t)/dt ^ 0, Eq. B5 is 
now approximated by 



dA^ 



Tf- 



(tr) 
FUV 



(t) 



dt 



-Ao 



e FUV 



(t)/Ao 



(B6) 



(Ao = (tNo), with the solution of the latter being 



Ag{jV(t) = Aoln [1 + (e<V(o)/Ao _ i^^-t/r, 



(B7) 



The latter expression reduces to that in Al when No — 00 (Ao 00), as expected. 
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Fig. 1. — The equilibrium molecular fraction fm. Plotted are contours (labelled with the corre- 
sponding value) of fm for various representative combinations of formation rate parameter fx, cloud 
power law index a, metallicity Z, radiation field Go and velocity dispersion Vdisp (Actual values 
are indicated in the plots). Above 1000 K the formation rate is taken to be zero, a indicates the 
power law index of the density profile (n oc r~"), thus a = means constant density clouds, a = 1 
logotropic clouds. 
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Fig. 2. — Comparison of the solution of Eq. 12 for the photodissociation of a molecular (plane- 
parallel) cloud and the solution of Goldshmidt &: Sternberg (1995) of the full integro-differential 
equations. Plotted is the total HI column N(HI) as a function of time since the onset of FUV 
irradiation (normalized on r/) where the drawn lines are the solutions of our Eq. 12 for different 
values of the Goldshmidt &; Sternberg (1995) a parameter (q;$ = rdis) while the dashed lines 
correspond to the solutions as read off from their Figure 1. Lines are labelled with the corresponding 
a parameter. Here we have adopted the same values for a a and $ as Goldshmidt h Sternberg (1995). 
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Fig. 3. — Overview of the ISM model for Z = Z0/5 (a typical metallicity for dwarf galaxies): 
equilibrium plots of (a) pressure P, (b) temperature T, (c) cooling rate A and (d) electron fraction 
Xe as a function of density n for three different values of the UV field Go (given in units of 1.6 x 10~^ 
erg cm~^ s~^) 
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Fig. 4. — Star formation of the model galaxies. Left panel shows the star formation as a function 
of time (Dotted line: high Z models, drawn line: low Z models). The right panel shows the 
mean star formation density as a function of radius (again for both models). The dotted lines in 
the right panels show exponentials with scale lengths of 0.7 and 0.75. 
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Fig. 5. — H2 fraction vs time. Dashed line: the molecular fraction fm for the ^ = 3.5, Z = Zq/5 
simulation, dash-dotted Une: fi = 17.5, Z = Z0/5, drawn line: fi = 3.5, Z = Zq, long-dashed 
line: /j, = 17.5, Z = Zq, and Dotted line: fi = 3.5, Z = Z0/5 with logotropic clouds. Simulations 
are started with fm = 0. 



-38- 





Fig. 6. — Radial dependence of the molecular fraction. Drawn line: mean molecular fraction as a 
function of radius (scale on left y-axes), and dotted line: total gas surface density (scale on right 
y-axes). Note that the y-axis scales are not the same for the different panels. 
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Fig. 7. — z-Height dependence of the molecular fraction. Drawn line: mean molecular fraction as 
a function of height above the disk plane (scale on left y-axes), and dotted line: total gas mass 
distribution (scale on right y-axes). Note that the y-axis scales are not the same for the different 
panels. 
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Fig. 8. — HI and H2 maps. Left panels: results for Z = Zq/S, right panels: Z = Zq. Shown 
are, from top to bottom, the HI distribution, the H2 distribution for = 3.5 and for ^ = 17.5. The 

HI maps of low and high ^ simulations arc very similar, shown are only the maps for /x = 3.5. The 
H2 maps have been divididcd by the global molecular fraction such that the differences with the HI 
distribution can be seen more clearly. 
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Fig. 9. — Power spectra of HI maps. Thick lines give 1-D power spectra of the LMC HI distribu- 
tion, averaged over all image lines in north-south or east-west directions (Elmegreen et al. 2001). 
Thin drawn lines are the corresponding power spectra made from projected HI maps of a simu- 
lated dwarf galaxy (Z = Zq, = 3.5). Given are the mean powerspecrum of 10 snapshots spanning 
200 Myr simulation time and the two lines indicating the mean ztla. Dotted lines are for the 
same model at low resolution (N=20k). The power is arbitrarily scaled, but the low and high 
resolution models are on the same scale and the LMC spectrum is scaled to match. 



Fig. 10. — Ha and H2: simulated Ha maps with H2 contours. Left panel shows the maps for the 
Z = Zq/5, /U = 3.5 simulation (contours at column densities of Nhi = 2 x 10^^ and lO^'^ cm~^) while 
the right panel shows the Z = Z©, /x = 17.5 run (contours at column densities of Nhi = 2 x 10^° 
and 10^^ cm~^). 
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Fig. 11. — H2 fraction vs time for simulations with different fgf. Dotted line: fgf = 2.5, drawn 
line: fgf = 10, dashed line: fgf = 20. The dash-dotted hne gives the molecular fraction for a 
star formation recipe using a threshold fm,sf = 1/8 (/x = 3.5 and Z = in all cases). 
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Fig. 12. — Distribution of star formation timescales for molecular regulated star formation. Shown 
are the cumulative distributions of star formation delay times, as measured from the time a particle 
becomes unstable and normalized on the local free-fall timescale tg, for different values of the 
treshold molecular fraction fm,sf- Dash-dotted line: im,si = 1/8) drawn line: ira,sf = 3/8 and 
dashed line: fin,sf = 5/8 (For comparison the dotted line gives an exponential distribution). 
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Fig. 13. — H2 fraction vs time for simulations with cooling. Drawn line: Z = Zq/5, dashed 
line:Z = Zq/5 with H2 cooling, dash-dotted hne: Z = Z© and dotted: Z = Z© with H2 cooling. 
For all: fi = 3.5. 
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Fig. 14. — Fraction of H2 vs total column densities. Gray scales give the distribution of H2 fraction 
NH2/NHI and total column density Nhi+H2 of the pixels of HI and H2 projection maps. Histograms 
give the distribution of pixel values of the total column density, scaled to the maximum bin. 
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Table 1: Overview of the processes included in the ISM model used. For H and He ionization 
equilibrium is explicitly calculated, for other elements collisional ionization equilibrium (CIE) is 
assumed, references: 1: Wolfire et al. (1995), 2: Raga et al. (1997), 3: Verner & Ferland (1996), 
4: Silva & Viegas (2001) 



process 


comment 


ref. 


heating 






Cosmic Ray 


ionization rate CcR = 1-8 10~^^ 


1 


Photo Electric 


FUV field from stars 


1 


cooling 






e,Ho impact 


H,He,C,N,0,Si,Ne,Fe 


2,4 


ionization 






& recombination 






UV 


ionization assumed for 






species with Ei < 13.6eV 




Cosmic Ray 


H, He only; primary 


1 




&: secondary ionizations 




Collisional 


H, He only 


3 


Radiative recombination 


H, He only 


3 


CIE 


assumed for metals 





